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Among the most popular and well studied quantum characterization, verification and validation techniques is 
randomized benchmarking (RB), an important statistical tool used to characterize the performance of physical 
logic operations useful in quantum information processing. In this work we provide a detailed mathematical 
treatment of the effect of temporal noise correlations on the outcomes of RB protocols. We provide a fully 
analytic framework capturing the accumulation of error in RB expressed in terms of a three-dimensional random 
walk in ’’Pauli space.” Using this framework we derive the probability density function describing RB outcomes 
(averaged over noise) for both Markovian and correlated errors, which we show is generally described by a 
gamma distribution with shape and scale parameters depending on the correlation structure. Long temporal 
correlations impart large nonvanishing variance and skew in the distribution towards high-fidelity outcomes - 
consistent with existing experimental data - highlighting potential finite-sampling pitfalls and the divergence 
of the mean RB outcome from worst-case errors in the presence of noise correlations. We use the Filter- 
transfer function formalism to reveal the underlying reason for these differences in terms of effective coherent 
averaging of correlated errors in certain random sequences. We conclude by commenting on the impact of these 
calculations on the utility of single-metric approaches to quantum characterization, verification, and validation. 


I. Introduction 

Quantum verification and validation protocols are a vital 
tool for characterizing quantum devices. These take many 
forms M, but one of the most popular due to its efficiency 
is randomized benchmarking (RB). In this approach, devel¬ 
oped originally by Knill et al. El, and expanded theoretically 
by various authors inoHia, the average error probability of 
a quantum gate {e.g. a bit flip) is estimated by implementing 
a randomly sampled gate sequence from the set of Clifford 
operations, and measuring the difference between the ideal 
transformation and the actual result. Averaging over many 
randomized sequences yields information about the underly¬ 
ing gate fidelity. 

RB has become so important in the experimental commu¬ 
nity Emma that despite the experimental complexity it is 
now common to simply quote a measured gate error, prb, re¬ 
sulting from a RB measurement for a particular experimental 
system, and relate this number to tight bounds such as fault- 
tolerance error thresholds in quantum error correction El- 
Of late, reported values of prb have compared favorably with 
these thresholds, and have been used to justify the scalability 
of particular experimental platforms. Underlying this entire 
approach is the subtle but central assumption that errors are 
uncorrelated El: an assumption which is generally violated 
by a wide variety of realistic error processes with long tempo¬ 
ral correlations iflTl l25l l26l . Violation of this assumption has 
been noted previously as a cause of increased variance of out¬ 
comes over randomizations, but until now there has not been 
a quantitative means to understand the impact of such tem¬ 
poral correlations, a detailed physical mechanism to explain 


why such distortion of RB results can appear, or a clear under¬ 
standing of the impact of this observation on the applicability 
ofRB. 

In this manuscript we examine the impact of relaxing the 
Markovian-error assumption by studying its effect on the dis¬ 
tribution of measurement outcomes over randomizations. We 
find that while all randomizations (meeting certain criteria) 
are valid within the RB framework, they exhibit starkly dif¬ 
ferent susceptibility to error when those errors exhibit tempo¬ 
ral correlations over multiple individual Clifford operations. 
We provide a detailed mathematical treatment of error accu¬ 
mulation in RB, using a general model treated in the specific 
case of dephasing errors. We demonstrate how the reduction 
of fidelity for a particular sequence may be given a geomet¬ 
ric interpretation, taking the form of a random walk in a 3D 
Cartesian coordinate system. The steps of this walk corre¬ 
spond to the appearance of Pauli X,Y ,Z errors derived from 
interleaved dephasing operators in the sequence of Clifford 
operations, and the overall statistics are determined by the un¬ 
derlying correlations in the noise process. Our treatment in¬ 
cludes both extremal cases of uncorrelated (Markovian) and 
systematic (DC) errors, as well as generic correlations inter¬ 
polating between these limiting cases and captured through a 
noise power spectral density. Our results provide simple an¬ 
alytic forms for the probability density functions of fidelity 
outcomes over RB randomizations as a gamma distribution, 
building connections to engineering literature on failure anal¬ 
ysis. We describe the impact of these observations on the 
interpretation of RB measurement outcomes in various noise 
environments and highlight the disconnect between measured 
RB error and metrics relevant to fault-tolerance such as worst- 
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case errors. 
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II. Randomized Benchmarking Procedure 

A common experimental implementation of single-qubit 
RB involves measuring the fidelity of net operations com¬ 
posed from random sequences of Clifford operations. For a 
sequence of length J, the net operation is written as the prod¬ 
uct 


i=i 


Vj 


( 1 ) 


of Clifford operators (see Table IIIi indexed by the se¬ 
quence T] = where me are random vari¬ 

ables uniformly sampled from the set {1,2,..., 24} labelling 
the elements of the Clifford group. Typically one makes the 
constraint on r/ that, absent any error processes, the ideal RB 
sequence performs the identity operation Sri = I. Of the 24*^ 
total sequence combinations for a given J, only a small sub¬ 
set of fc <C 24*^ random Clifford sequences is implemented 
in practice. Each sequence is repeated (with appropriate qubit 
reinitialization) n times to build an ensemble of fidelity mea¬ 
surements sampling and the underlying error process. These 
measurement outcomes are then averaged together to obtain 
useful information from projective qubit measurements, with 
the mean hdelity retained for each random sequence. 

We represent the space of measurement outcomes by the 
(fc X n) matrix 


rii / Fi^i Fi^2 




V2 


Vk 


f^2,l f^2,2 


.. Fi,„\ 

• ■ f^2.n 


( 2 ) 


\ Fk^l Fk^2 ■ ■ ■ Fk^n / 


where element Fij G [ 0 , 1 ] is the measured hdelity when im¬ 
plementing the ith Clifford sequence r/^ = (ryip, Vi,j) 

in the presence of the jth realization of the error process, de¬ 
noted by dj. For clarity, we make the distinction between the 
measured hdelity Fij (e.g. obtained in experiment via pro¬ 
jective qubit measurements) and the calculated hdelity (see 
Eq. 1^ serving as a proxy for Fij in our analytic framework. 
The process of repeating each rjj and averaging the hdelity 
outcomes over the finite sample of noise realizations therefore 
corresponds to a measurement of the “noise-averaged” hdelity 
for the jth sequence, represented by the quantity 


1 ^ 

{!,...,fcj 


( 3 ) 


i=i 


where the angle brackets (•) in the second subscript denote av¬ 
eraging over the columns of F^'^\ Similarly the mean hdelity 
averaged over both Clifford sequences and errors is denoted 






kn 


( 4 ) 


Z = 1 j = l 


where angle brackets in both subscripts indicate averaging 
over both rows and columns (that is, all elements) of F^‘^\ 
The quantity is therefore an estimator of the true mean 
hdelity = {F)n.s, formally obtained as the expectation 
over all possible hdelity outcomes F dehned on the support 
of the random variables rj and 6. 

In the standard RB procedure, measurements of fSF for 
increasing J are htted to an exponential decay from which 
the mean gate error prb is extracted as the decay constant. 
Implicit in this procedure, however, is the assumption that, 
for any random set of Clifford sequences, the resulting dis¬ 
tribution of hdelity outcomes represents the underlying error 
process fairly, and the total mean is reasonably repre¬ 
sentative of any individual sequence. That is, the distribution 

of values is symmetric about fi^F with small relative 

variance for any random set {rh ,..., 77 ^,}. 

In this paper we show to be an unbiased and effec¬ 
tive estimator only when the error process is truly Markovian. 
That is, it possesses no temporal correlations. We provide a 
detailed analytic calculation of the impact of realistic, corre¬ 
lated noise on the distribution of outcomes over different ran¬ 
domizations. To examine this we introduce a physical model 
compatible with the experimental procedure, and which per¬ 
mits accounting for the presence of noise correlations. 


III. Physical model 

In this section we develop the mathematical framework 
used to model and investigate the impact of noise correla¬ 
tions on RB. Our model assumes a Unitary error process with 
temporal correlations, rehecting typical experimental condi¬ 
tions where qubit rotations, e.g. from a fluctuating classical 
field, dominate. This process is generically described by a 
power spectral density (PSD) capturing the various correla¬ 
tion timescales. Our error model is fully general and considers 
universal (multi-axis) errors. For simplicity of technical pre¬ 
sentation we treat single-axis dephasing in the main text and 
provide the full presentation of universal noise in Appendix 

IS 

The main result of this section is the geometric interpre¬ 
tation of error accumulation in RB sequences in terms of a 
random walk in 3 dimensions (a step taken for each Clifford), 
with step lengths and directionality inheriting the correlation 
structure of the error process. This result is independent of 
any choice of representation of the Clifford group, though our 
numeric simulations later used to verify it do rely on a partic¬ 
ular representation (see Table [nl)i. 

A. Noise-averaged fidelity {F) 

We begin by introducing our metric to quantify the reduc¬ 
tion in fidelity for a given r/ and 6. We employ the trace h- 
delity 


F{rj,6) 



( 5 ) 


capturing the overlap between ideal, Sr], and noise-affected 
sequences, Sri,s, via the Hilbert-Schmidt inner product. For 
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inputs rjj and 5j this provides us with a computational metric 
corresponding to measurements Fij obtained in experiment. 
Consequently our proxy for the measured noise-averaged h- 

delity takes the computational form 


As we will show, the noise interaction “steers” the operator 
product away from the identity gate performed by the ideal 
sequence 5^ and reduces the operational hdelity in some way 
characteristic of the correlation structure of 6. 





C. Analytic Expression for Sequence Fidelity 

1. Series expansion for Sequence Fidelity 


where denotes an ensemble average over n realizations 
of 5, and explicit reference to t] and S in the argument of F 
has been dropped. Henceforth we also drop the subscripts 
and denote the calculated noise-averaged fidelity simply by 
(F). In the following sections we proceed with our main task; 
deriving the probability density function (PDF) of (F), and its 
dependence on the correlation structure of S. 

B. Error model 

To reiterate, a Clifford sequence r/ implements the opera¬ 
tion Srf = n/=i where C,,. is the p* Clifford generator. 
The random variables rji are uniformly sampled from the set 
{1,2, ...,24}, and are assumed to be independent and iden¬ 
tically distributed (i.i.d.), subject to the technical constraint 
Sri = I. Unitary errors are implemented by interleaving 
Sri with a sequence of stochastic qubit rotations, yielding the 
noise-affected operation 

Sn,s = UiCri^U2Cn^---UjCjij. (!) 

This approach holds for arbitrary unitaries, Uj, enacting rota¬ 
tions in any Cartesian direction; a full treatment for universal 
noise models appears in Appendix [D| However, in the follow¬ 
ing we restrict ourselves to the subset of unitaries enacting 
a sequence of dephasing rotations, Uj = eytp{—i5jZ), pa¬ 
rameterized by the time-series vector 6 = { 61 , 62 , ■■■, 6 j). We 
assume the error process is wide-sense stationary {i.e. the for 
errors 6 ^ and 6 j the two-point correlation function { 6 i 6 j) = 
{6i6i-k) depends only on the time difference k = i — j), and 
has zero mean. 

Temporal noise correlations are established by introducing 
correlations between the elements of 6. In this work we treat 
three distinct cases 

(a) Markovian process: Elements 6j ~ A/'(0, cr^) are i.i.d. 
Gaussian-distributed errors with zero mean and vari¬ 
ance variance cr^, and completely uncorrelated between 
distinct Clifford gates in any sequence Sri (correlation 
length 1). 

(b) DC process: Elements 6j = 6 ^ A/'(0, cr^) are identical 
over a given sequence 5^ (maximally correlated, corre¬ 
lation length J), but are i.i.d. (uncorrelated) Gaussian- 
distributed errors with zero mean and variance variance 

over different instances. 

(c) Generically-correlated process: Correlations between 
elements of S separated by a time interval of “k gates” 
in Sri are generically specified by an autocorrelation 
function Cs{k) = { 6 j 6 j+k), in terms of which S may 
be described by a power spectral density (PSD) S{uj). 


We begin by obtaining an approximation for {F). As RB is 
intended to estimate errors too small to resolve in a single gate 
implementation, we assume the noise strength per gate, cr, is 
small. Over long gate sequences the accumulation of errors 
may be quantified by Ja^, which can be large for sufficiently 
large J. Assuming Ja^ < 1, however, we can make analytical 
progress by expressing the error unitaries in terms of a power 
series truncated at O (cr^). The noise-affected sequence Eq. 
is therefore approximated by 


« n j ( 8 ) 

The full expansion of Eq. generates products of order up 
to o(ni=i^i) = Let denote the 

O ^ (^") product due to cross-multiplying 

terms like {6j^Z)^^ ...{6j^Z)^’'^, where = 

n. Retaining terms only up to fourth order (n = 4), consis¬ 
tent with the order of our original Taylor approximation, we 
thereby obtain 


5, 


Tj.S 


«‘”'+4'>+ca 


The first term in this expansion is identical to the ideal Clif¬ 
ford sequence, Higher-order terms capture suc¬ 

cessive error contributions and are composed of blocks, or 
subsequences, within Sri interrupted by one or more Z er¬ 
rors, indicated by the subscripts. To evaluate Eq. we 
must obtain expressions for the quantities = 

iTr k )■ Detailed derivations for these terms are 

given in Appendix [C] 


2. Mapping from Clifford Sequence to “Pauli space” 

We compute the terms in the power series by relating a 
given Clifford sequence to a random walk in “Pauli space”. 
Eor illustrative purposes, we consider in detail only the 
quadratic terms: 

= - II 6j6kCi,j-iZC,,k-iZCk,j (10) 

j<k 

where the ordered summand runs over 1 < j < k < J, and 
we have dehned the Clifford subsequence operators 


Cjk = Cn^-Cn„ 1<J <k<J. (11) 
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We now define the cumulative operators Km giving the prod¬ 
uct of the first m Clifford operations 


Km — Cim — 


( 12 ) 


where Kq = Kj = I, and = Km^ is also a Clifford 
operator (since it is a product of Clifford operators). Any sub¬ 
sequence Cjk therefore “factorizes” as 


Consequently, i? is a random walk in “Pauli space”. 

Taking half the trace of Eq. and substituting in Eq. 20 


the term cancels out in the expression for ^Tr ( Sri^s ) ■ 
This further simplifies by observing that 


1 . 


■'( 




0, k odd 
1, k even 


( 21 ) 


Cjk — (13) 

allowing us to rewrite 

= P,Pfe (14) 


This follows from (a) the cyclic property of the trace, (b) the 
technical constraint Srj = Cij = 1, and (c) is either Z or 
I depending on whether k is odd or even respectively. Using 
Eq. [ 21 1 we hnd = 1 and = 0, 

leading to the simplihed expression 


where the operators on the right hand side are dehned by 

Fm = Km-iZKl_, G {±X, ±Y, ±Z} (15) 

for 0 < m < J, and are always signed Pauli operators (since 
the Clifford group is the normalizer of the Pauli group). We 
may therefore express the operators P^ in the basis of Pauli 
operators as 


Pm — VniF Zm^ 


(16) 



«l-i|lH|l2 + Q® ^ (22) 

+ + Qhl,2 + Q^ll + * 32^2 + * 34 "^^- 


Substituting this into Eq. and retaining only terms up to 
O we obtain 


(J-)«l-(||^f)+0(4)^ (23) 


where 


where Xm,ym,Zm G {Ojil} subject to the constraint 
\xm\ + \ym\ + \zm\ = 1 (only One nonzero coefficient). 
The unit vector defined by 

^m = (^mj ym; ^m); Ijl’mll — 1 (1^) 


therefore points uniformly at random along one of the princi¬ 
ple Cartesian axes {±a;, ±y, ± 2 :}, and maps the “direction” 
of the operator Pm in Pauli space. Thus we have constructed 
a map from a given random Clifford sequence of length J, 
to a set of J unit vectors each oriented at random along the 
cartesian axes. 


With these insights the error contributions may be recast 
into more convenient expressions by moving to vector nota¬ 
tion. In particular, taking the trace over Eq. 14 and using the 
cyclic composition properties of the Pauli matrices, we find 


^Tr(PjPfc) = • ffc. (18) 

Prom Eqs. [^[^and[T^we therefore obtain 

Q^ll= ■ Tk- (19) 

j<k 


ow = i (^\\R\\^) + 2 + 2 (gf ) . (24) 

To a good approximation may be treated as a small cor¬ 
rection in the form of a small constant, with the statistical 
distribution properties residing in the leading-order term, the 
random variable (||i?|p). In arriving at these expressions we 
have used the assumption that the error process has zero mean, 
from which it follows only terms with noise random variables 
raised to even powers, or those summed over terms raised 
to even powers, survive the ensemble average with the oth¬ 
ers reducing to zero. In particular, (Q[^i i) = 11 ) = 

(* 3 m. 2 ) = (* 3 i'^ 3 ) = 0 (see Appendix [c|. 

IV. Probability Distribution Functions for RB Outcomes via a 
Geometric Interpretation of Error Accumulation 

In the expressions above we see that the key metric cap¬ 
turing the reduction of fidelity in a randomized benchmarking 
sequence is 


J 




(25) 


Observing the quantity SjSkVj ■ Vk is invariant under ex¬ 
change of indices, we may recast the restricted sum over 
j < k into an unrestricted sum, picking up a residual term 
= — I J2j=i (see Appendixj^, to obtain 

= R=Y.^3r,. (20) 


which may be interpreted as a random-walk in “Pauli 
space” generated by adding J randomly-oriented steps along 
the principle Cartesian axes, with step lengths specihed by 
8. Walks which terminate far from the origin correspond to 
sequences with large net infidelities while those which ulti¬ 
mately end near the origin have small infidelities. These ob¬ 
servations form a key contribution of this work, as we will 
show. 
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This geometric picture provides a unique insight into how 
error accumulates in RB sequences, and facilitates calcula¬ 
tion of the distribution of {T) incorporating the effects of cor¬ 
relations in the error process. In particular, the calculation 
of sequence fidelities maps onto the distribution of random 
walk, with possible correlations manifesting the random step 
lengths. 

In the following sections we explore how correlations in 
S affect the distributions of the terms in Eq. and hence 
the probability density function (PDF) (F) of the noise- 
averaged fidelity (F). Our calculations demonstrate that un¬ 
der all error models treated here (F) takes the form of 
a gamma distribution, which is well known in statistics and 
provides the significant benefit in that explicit analytic forms 
are available for both the moments of the distribution and the 
moment-generating functions (see Table |I]). Interestingly the 
gamma distribution is known to be useful for failure analysis 
and life-testing in engineering lIZTl which bears some similar¬ 
ity to the notion of error accumulation in RB. 

A. PDF for Markovian Processes 

In the Markovian limit (i.e. uncorrelated noise), we assume 
all noise random variables dj are i.i.d. Hence R corresponds 
to a J-length unbiased random walk with step lengths sam¬ 
pled from the normal distribution A/(O, cr^). Since these 
step lengths are symmetrically distributed about zero, the 
distributions of the components of the walk vector SjVj = 
{SjXj, SjUj, 6jZj) are invariant with respect to the sign of the 
coefficients aj in all Cartesian directions a G {x,y,z}. Ig¬ 
noring the signs we therefore treat the the coefficients as bi¬ 
naries aj G {0,1}, where the zero event simply reduces the 
number of steps taken in that direction. Let 


,7 

n„ = ^|aj|, aG{x,i/,z}, + Uy + = J (26) 

i=l 

count the total number of nonzero components in each 
Cartesian direction over the sequence of walk vectors 
...,fj}. Then 

R=(^s^ + ... + s:^, sf + ... + sy^, S^ + ... + Sl) (27) 

where the superscripts in indicate summing only over the 
subset of Sj for which the coefficients aj are nonzero. Thus 
we have 

\\Rr = Al + Al + Al (28) 

where we have defined -f + ■•■ + for 

a G {x,y,z}. Since all Sj ~ Af (O, a^) are i.i.d., the new 
random variables Aq, ^ M (O, UaCr^) are also independent 
and normally distributed, with zero mean and variance nacr'^. 

The distribution of the sum of squares of non-identical 
Gaussians is generally complicated to write down, requiring 
a generalized chi-square distribution. To avoid this we make 
the following modest approximation. Since the vectors Vj are 
uniformly-distributed there is a f probability of being parallel 


to any given Cartesian axis. The probability of finding any 
particular combination {nx,ny, Uz) is therefore given by the 
multinomial distribution 


V {nx,ny,nz) 



(29) 


For J > 5, however, this is sufficiently peaked around 
'n-x,y,z = J/i that we may simply regard these values as fixed 
without significant error. In this case A^^y^z Af (O, jS) 
reduce to i.i.d. random variables. The distribution of ||ii|p 
consequently reduces to chi-square distribution with 3 degrees 
of freedom. It is more convenient, however, to express this in 
more general terms as a member of the two-parameter family 
of gamma distributions (see Eq. |A5[ ), of which the chi-square 
is a special case. Specifically, we obtain 


ll^f ~r(a,/3), a=^, 

with shape parameter a and scale parameter ji. An ensem¬ 
ble average over n independent noise realizations is therefore 
specified by 


\Rf)n = -J2\\R\\l ||B 


i=i 


2’ 3 


(31) 


where the ||.R||| are i.i.d. gamma-distributed random vari¬ 
ables. But the sample mean over n gamma-distributed ran¬ 
dom variables simply yields a rescaled gamma distribution 
with a —?■ na and f] ^ (3/n (see Eq. A7 1 . Consequently 


{\\Rr)n- 

with expectation E[(||.R|p)„] 
V[(|l.Rf )„] = 


3n 2Ja^ 
T’ 3n 


(32) 


= JtT^ and variance 


Higher-order contributions may be included by computing 
the terms in For full derivations see Appendix]^ here 
we sketch the result. From the known distribution of |j.R|p 
the PDF for ||i2||‘^ is given by a transformation al- 


E6 


in Eq 

lowing us to compute the expectation and variance. Taking an 
ensemble average over noise realizations, applying the cen¬ 
tral limit theorem, and observing the narrowness of the re¬ 
sulting distribution compared to the leading-order term, we 
make the approximation (||.R||'*) « E[||.R||^] = The 

remaining terms yield values (< 32 ^ 2 ) = ~ 

we obtain the 


Substituting these intoEq. 


4th order correction ~ P 
fidelity reduces to 


2^4 


24 


and the noise-averaged 


{T)^l-{\\R\\P + -J^a^ 


(33) 


inheriting the gamma distribution of 
Eq. p 


i2||2)„ described in 
We linearly transform this expression to produce a 
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final PDF for noise-averagedin the Markovian regime: 

f^r) (F) = (34) 

where r(a;)is the gamma function, and the quantities I'iF), a 
and P are defined in Table |I] 

B. PDF for DC (quasi-static) processes 

In the DC limit (i.e. quasi-static noise) we assume all noise 
random variables 6j = (5 are identical (maximally correlated) 
over a given sequence Srj- However over separate instances S 
is sampled from the normal distribution S ^ JV {O, . Then 

R corresponds to a J-step unbiased random walk with fixed 
step length 6 directed along the Cartesian axes. In this case 
the noise random variables 6 and Clifford-dependent random 
variables rj factorize, allowing us to express 



Markovian 

DC 

B lock-correlated 

a 


3 

2 

fj/(M-l) 

0 

/n 

|JF 


u{F) 

1 - F-H 

1 - F 

1 - F 

E 

1 — JF + fFF 

1 — J(T^ 

1 — 

M 

I-JF (l-^) + |jV4 

\-\jf 


V 

jn 


|J(M- 1)0-4 

S 

-2V2/3n 

-‘i-s/Ffii 

-2^2(M - 1)/3J 


TABLE I. Scale and shape parameters for the respective gamma dis¬ 
tributions, and calculated moments for noise-averaged fidelity dis¬ 
tributions (F). Moments obtained from standard gamma dis¬ 
tribution after appropriate linear transformation specified by viF): 
Expectation, E = u + F — aj3'. Mode, M = u + F — 
{a — l)/3; Variance, V Skew, S 


R = SV, V = Y^ Tj (35) 

i=i 


where V G defines an unbiased random walk on a 3D lat¬ 
tice generated by adding J unit-length steps. Thus, in contrast 
with the Markovian case, here the random walk in Pauli space 
is unaffected (step-by-step) by the noise interactions. Rather, 
in a given run, the noise variables S effectively scale the ran¬ 
dom walk V generated by the Clifford sequence, up to a sign. 
Since we are interested in the norm square ||i?|p = V|p, 
however, any sign dependence of S vanishes. Performing a fi¬ 
nite ensemble average over n noise randomizations we there¬ 
fore obtain (||.R|p)„ = (^^)n||^|l^, and Eq. [j^yields 

(J^)«l-(52)„||Ff+ 0(4) (36) 


In this case, 0(4), includes a term j(^^)n||^||^ which is now 
highly correlated with the leading-order contribution, so we 
cannot use the expectation value as a proxy for the whole dis¬ 
tribution. However corrections from these terms are so 

we ignore these terms, and formally study the limit Ja^ +; 1. 
Numerical evidence indicates that this approximation works 
well up to Jcr^ ^ 1. 

Since <5 ^ Af (O, cr^) is normally distributed, 6'^ is chi- 
square-distributed which is a special case of the gamma distri¬ 
bution, (5^ ^ r(a, /?) with shape parameter a = 1/2 and scale 
parameter = 2a^. Taking the sample mean over n inde¬ 
pendent gamma-distributed variables again yields a rescaled 
gamma distribution with a —>■ na and ^ fi/n 



with expectation E[(i5^)ri] = and variance = 

2+) 

n 


The random-walk behavior of || V|p (Eq. 351 represents a 
well-studied problem in diffusion statistics. Let the random 
variable R be the distance from the origin in a symmetric 
(Bernoulli) 3D random walk after J steps. It is straightfor¬ 


ward to show that the PDF for R is 

3/2 


= (sif) 


9 -3r^ 

47rr e . 


(38) 


This expression describes a random walk of unit step length, 
and is derived assuming that Vj is uniformly and continuously 

is a good 
The 


38 


sampled from all directions in however Eq. 
approximation for the PDF of a walk on a 3D lattice 
distribution of the distance square || « R^ is then given 


by the transformation (see Eq. A3 1 




2a;i/2- 

1 

r(a)/3‘ 


-x“ 4 exp 


(39) 

(40) 


where a = 3(2 and = 2J/3, and r(a:) is the gamma func¬ 
tion; again this is a gamma distribution (see Eq. |A5| l with 
shape parameter a and scale parameter /3. Consequently 


IVII 


r,'3 21 

2’ 3 


(41) 


Thus, to leading order, the PDF for {F) with DC noise 
is specified by the product of two independent gamma- 
distributed random variables. The closed-form expression can 
be calculated by direct integration (see Appendix]^, however 
for moderate ensemble sizes n > 50 it is sufficient to approx¬ 
imate {5'^)n as strongly peaked around its mean, tr^, such that 
{F) « 1 —CT^||V|p. In this case the PDF(F) reduces to a 

linear-transformed gamma distribution associated with || 
and possesses the same form as Eq. 34 but with different val¬ 
ues of the parameters 12 (F), a and /3, as given defined in Ta¬ 
ble 1^ This is a remarkable observation given the substantial 
differences in error correlations between these extreme limits. 

C. PDF for Generically-Correlated Processes 

For both limiting cases of Markovian and DC noise cor¬ 
relations treated above we showed (F) is described, to first 
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order, by the gamma distribution with shape and scale pa¬ 
rameters dependent on the correlation structure. Assuming 
continuity of the distribution we also expect (T) to be ap¬ 
proximately gamma-distributed for generic, intermediate cor¬ 
relation structures. We can show this formally for a spe¬ 
cific class of correlated noise models, in which the noise 
is block-correlated', i.e. the error random variables 5j are 
constant over blocks, or subsequences, of Clifford operators 
of hxed length M < J, and there is no correlation be¬ 
tween distinct blocks. The full derivation of the PDF for this 
case can be found in Appendix We find block-correlated 
noise yields a gamma-distributed fidelity, with parameters 
a = ZJI2(M — 1), /3 = 2{M — l)tT^/3, enabling us to in¬ 
terpolate between the Markovian (M = 1) and DC (M = J) 
limits for arbitrary correlation length M (see Table|I]). 

While block-correlated noise is not stationary (i.e. it does 
not have a stationary power spectmm), it simply and explicitly 
captures the notion of a correlation length, M. The correla¬ 
tion length thus manifests itself in the distribution of fidelities. 
With these insights, for brevity we assume that generic noise 
correlations also give rise to gamma-distributed infidelities. 
This is supported by the quantitative calculations and qualita¬ 
tive arguments in Appendix [G| and by comparison with Monte 
Carlo numerics. Consequently, the shape and scale factors can 
be inferred from the mean and variance of the distribution, 
which we calculate directly. 

For simplicity, we assume the error process is sufficiently 
non-Markovian that the noise ensemble size n and the O (u^) 
contributions may be ignored without introducing large errors, 
as in the DC case. Thus, we write 

(||,Rf)^r[a,/3]. (42) 


From the hrst two moments of the gamma distribution, the 
expectation and variance of the random-walk variable are 
E[(||.R|p)] = a/3 and V[(||.R|p)] = a/3^, so the shape and 
scale parameters are given by 


a = 





(43) 


The task of defining the gamma distribution therefore reduces 
to computing the expectation and variance of (||i2|p), as a 
function of the generic correlation structure of the error pro¬ 
cess, S. 

An arbitrary time-series x(t) may be characterized by 
an autocorrelation function Cx(t) = {x{t)x{t -I- r))( 
where the ensemble average is over all t, and r is the 
time difference between measurements. In this case, in¬ 
voking the Wiener-Khintchine theorem, the power spectral- 
density (PSD) is given by the Fourier transform S{uj) = 
Cx{T)e^‘^'^ dr. We make use of these relations by 
discretizing the time-series. Let the elements of d be defined 
by 5j — x{tj) for tj/xg G {1, 2,.., J} where Xg is the time 
taken to perform a Clifford operation. The underlying error 
process is thereby discretely “sampled” by the Clifford se¬ 
quence Sr), and correlations between elements of S separated 
by a time interval of “k gates” are specified by the (discrete) 


autocorrelation function 


Cs{k) = {SjSj+k)j (44) 


where the ensemble average is over the index j. Substituting 
these definitions into Eq. [^the expectation and variance of 
(|j.R||^) are given by (see Appendix ^ 


. J-i 

E = JCs{0), Y = -J2(J-k){Cs{k)f. (45) 


In experiment one typically has access to the PSD S{uj) rather 
than the autocorrelation function. However this easily maps to 
our framework via an inverse Fourier transform 





(46) 


The PDF therefore has the same form as Eq. 34 but with a 
and P given by Eq. 43 and v{F) = \ — F. 

With these expressions we now have a complete analytic 
representation of the distribution of measured noise-averaged 
fidelities over different RB sequences. In the next section we 
verify these results with numeric Monte Carlo simulations and 
discuss the differences in distribution characteristics depend¬ 
ing on noise correlations. 


D. Comparing PDFs for Various Noise Correlations 

The analytic forms for (F) derived above serve as a 
tool to analyze the impact of temporal noise correlations on 
RB experiments. In all cases (Markovian, DC, and Intermedi¬ 
ate) the PDE is r(a, /3) -distributed, with differences captured 
in the values of the shape parameter a and scale parameter /3. 
This result is derived from statistics of a 3D random walk. 

Plotting the PDE (F) in Eig. immediately reveals 
substantial differences in the distribution of outcomes for the 
two limiting cases. While the distributions in both Marko¬ 
vian and DC cases yield the same value for the mean (the 
statistic currently used in RB protocols), the higher order mo¬ 
ments diverge significantly. Eor DC noise (F) is skewed 
towards high-hdelities and possesses a variance signihcantly 
larger than that for Markovian errors of equivalent strength, 
parameterized by the value of cr^. Averaging over a large 
noise ensemble results in the mode converging to the mean 
in the Markovian regime, but maintaining a fixed higher value 
of fidelity for the DC regime. By comparison, the variance 
and skew for Markovian noise diminish with increasing noise 
averaging, but remain fixed and nonzero in the DC case. 

To compare our analytic PDEs with the true distributions 
we directly simulate the fidelity outcomes associated with the 
metric in Eq. Eor a given J, we generate an ensem¬ 
ble {rii, of k random Clifford sequences. The hrst 

J — 1 elements in each sequence r|^ are uniformly and in¬ 
dependently sampled from the set {1,..., 24}, with the hnal 
element chosen such that the total operator product performs 
the identity Sr). = I. Eor each sequence rf^ we then generate 
an ensemble ..., di n] of n random and independent re¬ 
alizations of the error process, where each Sij is a sequence 
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FIG. 1. Analytic PDFs (J^) and simulated distributions for different 
noise regimes using prb ~ 2 x 10“"^. a, b) (F) calculated as a 
function of J for Markovian and DC regimes, normalized to unity at 
the mode for each J for clarity. White line indicates analytic mean, 
E[(J^)]; red dashed line indicates analytic mode, M[(J^)]. c, d) Com¬ 
pare analytics (solid tines) to numerically simulated histograms of 
distributions for various J, using n — 50 and no free parameters. 
Curves vertically offset by for clarity by multiples of 25 units. In¬ 
set c) Markovian distributions varying n=10, 50, 250 (black, blue, 
pink), corresponding to distribution in dotted box. e) Analytic PDF 
(F) (solid red line) and simulated distributions for quasi-l/cu noise 
regime with J = 200 and n = 3000. Comparative PDF for DC 
(left axis) and Markovian noise (right axis) shown as black dashed 
and solid lines respectively. Noise parameters chosen such that the 
mean error is 0.05. PSD shown in inset, constructed using Fourier 
synthesis as described in (29). 


of J noise random variables generated by Monte Carlo sam¬ 
pling from the appropriate correlated-error model. For Marko¬ 
vian and DC processes, this is fully described in Section [nl| 
For the intermediate case, random sequences d with a desired 
PSD S{uj) and autocorrelation function Cs{k) may be gener¬ 
ated by uniformly-phase-randomized Fourier synthesis as de¬ 
scribed in ll29) (see Appendix|I]i. 

For each pair rj^ and 6 ij the operator product in Eq. |7]is 
computed, using Table to determine the (2 x 2) unitary 
matrix representing each Clifford operator. For each pair the 
trace fidelity Fij = F{ri^, Sij) is then calculated using Eq. 
generating the array shown in Eq. simulating the measured 
fidelity outcomes. Averaging over columns, as in Eq. the 


array reduces to a column vector containing k noise-averaged 

fidelities one for each t]^, which we finally plot as a 

normalized histogram and compare against our analytic PDEs. 


Pig. [TJc) and (d) compares the numerically generated his¬ 
tograms of hdelity against the analytic results, Eq. 34 for 
Markovian and DC noise respectively. These are are in ex¬ 
cellent agreement for the different error processes considered. 
The characteristic long-tailed distribution peaked near high h- 
delities in the DC limit reproduces key features observed in 
recent experiments ini. Agreement with analytics is good 
(with no free parameters) for Ja^ < 1 , beyond which higher 
order error terms contribute to the distribution. 


Por correlated noise, similarly good agreement is obtained. 
We have validated the block-correlated noise model (not 
shown). For quasi-l/w noise, the fidelities are also well de¬ 
scribed by the gamma distribution, accounting for correla¬ 
tion length, as shown in Fig.[TJe). The dominance of low- 
frequency components in this PSD (see inset) skews (F) 
toward higher fidelities than the mean, but the presence of 
higher-frequency components (shorter correlation times) re¬ 
duces the variance and partially restores symmetry. In the 
limit of power spectra containing substantial high-frequency 
noise, small shifts due to higher-order terms in Eq. [^become 
important. 

The underlying physical reason for the differences in the 
two limiting error models is revealed by examining the 
filter-transfer functions G(w; J 7 ) for the various Clifford ran¬ 
domizations OOl . The noise-averaged hdelity is given by 
{F)oo « (1 — e~^ -^ 0 °° Giuj;ri)S{ui)u) quantifying the 

susceptibility to error over a given frequency band, and has 
demonstrated experimental applications ED. Using tech¬ 
niques outlined in Refs. |[30] [32 we calculate G(a;;rj) for 
10^ random sequences 77 , shown in Pig. In the low- 
frequency regime we observe variations over several orders- 
of-magnitude in the vertical offset of G(a;; 77 ) and also varia¬ 
tions in slope at higher frequencies, indicative of partial error 
cancellation and hence substantial variations in susceptibility 
to correlated errors. The corresponding distribution of hdeli- 
ties agrees well with /^jr) (F) for DC noise (see inset). 

These observations arise from the fact that some RB se¬ 
quences contain coherent, error-suppressing subsequences. In 
fact RB sequences bear resemblance ll^ to randomized dy¬ 
namical decoupling protocols known to suppress errors in cer¬ 
tain limits of correlated noise ES- This would lead to “ar- 
tihcially” small measured error for correlated noise, thereby 
increasing the variance and skew in (F) towards high h- 
delities for sufficiently low-frequency-dominated noise. Pur- 
thermore, this explicit link between the form of (F) and 
underlying symmetries in suggests we may downselect RB 
sequences using, e.g. the hlter function or appropriate entropy 
measures on the Clifford sequences to ensure coherent aver¬ 
aging properties are minimized. To the best of our knowledge 
this is the first quantitative validation of the mechanism under¬ 
lying the shifts in (F) for temporally-correlated errors, 
and the hrst application of the hlter-transfer function formal¬ 
ism for quantum control at the algorithmic level. 
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FIG. 2. Filter functions {G(a;; rj)} for ensemble of RB sequences. 
Dimensionless angular frequency lj normalized to the duration of a 
bit-flip operation, low-frequency (w < 1) noise susceptibility is 
captured by vertical offsets and slopes. Inset: Histogram of calcu¬ 
lated fidelities for {G(a;; t])} using S{uj) oc S{lj — 47r x 
for quasi-DC noise (scaled to correspond to a = 0.015) and overlaid 
with/(^)(F’). 


V. Discussion and Conclusion 

Our primary observation is that the form of {F) - and 
in particular the moments of the distribution - can exhibit 
strikingly different behaviors in different error regimes. This 
is true despite the fact that the expectation E[(J^)] = 
approximately converges for Markovian and DC cases for 
weak noise (see Table |I]l. This has a variety of important im¬ 
pacts in the application and interpretation of RB experiments 
for quantum information. 

For instance, due to the form of /^jr) (F) for low frequency 
noise, the sample estimator obtained from an ensem¬ 
ble of A: ^ 24“^ instances of rj can differ from the true ex¬ 
pectation E[J^]j, generally leading to overestimation of the 
mean fidelity ifT^ fTSll . This risks systematically underesti¬ 
mating prb due to insufficient sampling over rj. The differ¬ 
ence — E[J^] j I may be formally bounded using moment¬ 
generating functions for the gamma distribution (see Supple¬ 
mental Material) to provide a more inclusive bound than pre¬ 
vious approximations imiia. Ensuring falls within an 
acceptable confidence interval, say within ±10% of E[±']j 
(and assuming a = 0.015, or ^ 2 x 10“^) requires 
> 9 randomizations be selected in the Markovian case, 
but > 443 for DC case. These values increase rapidly 

as confidence bounds tighten. These sample sizes are much 
larger than typically employed in experimental settings, but 
may be partially relaxed when calculating prb by fitting to 
measurements performed for multiple values of J. In the 
Markovian regime we also find a tradeoff between finite noise 


sampling (experimentally achieved by repeating a fixed se¬ 
quence Sr^ and averaging over the resulting projective mea¬ 
surements), and the required . Increasing the noise aver¬ 
aging, n, reduces the value k^^'^. 

Beyond the question of how well measurements performed 
in strongly correlated environments estimate prb, there is un¬ 
certainty surrounding the breadth of applicability of this single 
proxy metric as a general quantum verification tool for quan¬ 
tum information. One key observation is that while the vari¬ 
ance and skew of the (F) converge to zero for Markovian 
errors in the limit of infinite noise averaging (n —>■ oo), both 
remain fixed for DC noise. Accordingly our results provide 
direct evidence of the divergence between prb and parameters 
relevant to fault-tolerance 051 such as worst-case errors O^ 
in noise environments with strong temporal correlations; in 
the DC limit the worst-case error can be much larger than the 
average error. 

Finally, the fact that some RB sequences are intrinsically 
“blind” to correlated errors, highlights potential shortcomings 
in performing experimental gate optimization by maximizing 
p^F at fixed J; an operator may optimize an experimental 
parameter (to maximize the RB fidelity) in such a way to in¬ 
crease systematic errors in individual gates. Such issues may 
be partially mitigated by selecting subsets of valid RB se¬ 
quences using the length of the random walk, | |i?| |, in the DC 
limit to ensure fidelities are not overestimated in a RB proce¬ 
dure. Future work will explore the use of entropic measures to 
associate RB sequence structure with ||i?|| without the need 
to perform full calculations of the random walk. 

We conclude that the interpretation and applicability of a 
measured RB outcome prb can differ significantly depending 
on the nature of the underlying errors and measurement pa¬ 
rameters experienced in a real experimental situation. This 
challenge can be partially mitigated through presentation of 
more complete datasets - specifically F'i,(.) for each sequence 
- in order to assist readers making comparisons between re¬ 
ported results. We believe the new insights our calculations 
have revealed will help bound the utility of prb in quantum in¬ 
formation settings and also help experimentalists ensure that 
measurements are not subject to hidden biases. 
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A. Mathematical Preliminaries 

1. Linear transformation 
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Let y is a continuous, non-negative random variable described by probability density functions /y(j/). If Z is the random 
variable defined by Z = aV + /3, where a, /3 € K and a ^ 0, then the probability density functions of Z is 

fz(z) = r-[fy ( -(Al) 

|a| \ a J 

2. Product distribution 

Let X and Y be to independent, continuous random variables described by probability density functions fx (x) and fy {y)- Then 
the joint distribution of Z = XY is 


fz{z) 



1 


fx{x)fY{z/x)dx 


(A2) 


3. Strictly increasing fnnction of a random variable 

Let X be an absolutely continuous non-negative random variable described by probability density function fx{x). Dehne the 
transformed random variable Y = g{X) where g{x) is a strictly increasing function. That is Xi > X 2 g(xi) > g(x 2 ), so 

that g(x) has a well-dehned inverse g~^(x). Then the probability density functions of Y is 

fviy) = fxig~^{y)) ^^ , (A3) 

ay 

4. Sum of continuous random variables 

Let X and Y be two independent, continuous random variables described by probability density functions/x(a;) and fyiy). 
Then the probability density functions of the sum Z = X + Y is 

/ OO pOO 

fx{z - y)fY{y)dy = / fy^z- x)fxix)dx (A4) 

-OO J —00 


5. Gamma distribution 

The gamma distribution describes a family of continuous probability distributions, related to beta distribution, and arising nat¬ 
urally in processes for which the waiting times between Poisson-distributed events are relevant. The common exponential 
distribution and chi-squared distribution are special cases. The gamma distribution is a two-parameter distribution, and there 
are a few different parametrizations in common use. Throughout this paper we parametrize the distribution in terms of its shape 
parameter a and a scale parameter /3. Let ^ L (a, /3) be a gamma-distributed random variable, then the probability density 
function is defined by 


fx{x) 




Ot — 1 

-X exp 



where r(a:) is the gamma function. Table A 5 contains some relevant statistics. 


(A5) 


mean 

aj3 

variance 


skew 

2 /v^ 

mode 

{a - l)/3 


TABLE 11. Statistics of the gamma distribution r(a, /3). 


6. Distribution of sample mean of gamma-distributed random variables 

Let Xi ^ T{a,P), i G {1, ...,n}, be a set of n independent and identically distributed random variables sampled from the 
gamma distribution with shape parameter a and scale parameter /3. Then the sum Z = follows the transformed 

gamma distribution 


Z ~ r(na, /3) 


(A6) 
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The ensemble average Zjn has probability density function given by transforming nfzijix). Substituting into Eq. A5 
defining rescaled parameters a na and f3 —> /3/n, the distribution of the ensemble average {Xi)n is found to be 


and 


{X,)r 




(A7) 


7. Lindeberg-Levy central limit theorem 

Let {Xi,X 2 ,Xn} be a set of i.i.d. random variables with expectation E[2fi] = p, and variance V[2fi] = < oo. Then the 

sample average Sn = ^ converges to the normal distribution 

lim ^/n{Sn — p) ~ A/” (O, p^) (A8) 

n—)-oo ^ ' 



13 


B. Clifford Group Representation for Single Qubit 

A unitary operation C is an element of the Clifford group if 

= r 


(Bl) 


where we have defined the Pauli group 


V = {±l,±X,±Y,±Z}. (B2) 

That is, the Clifford group is the normalizer of the Pauli group, where for every Pauli operation P G V, there is another P' gV 
such that CPC^ = P'. For a single qubit, the set of all such C may be thought of as rotations of the Bloch sphere that permute 
the orientation of ±X, ±Y, ±Z in the Cartesian basis associated with the Pauli matrices, which we refer to as “Pauli space.” 
To obtain a clearer physical picture of these operations, consider associating X to any of the six Cartesian axes {±a:, ±y, ±z}. 
With this axis fixed, we may rotate the axes about X into four possible orientations while preserving XYZ right-handedness. 
This is the action of the C group; the group of rotational symmetries of the cube. 

We construct our representation as follows. Let Ri{d) represent one of nine elementary unitaries generating a clockwise 
rotation (looking down the axis of rotation toward the origin) through angle 0 G {tt, ±7 r/2} about axis i G {a:, y, z}. The three 
TT rotations correspond to 


Rx.y,z{Tr) = X,Y,Z 


(B3) 


and we use the shorthand 


Rf = P*(±7r/2), i G {x, y, z} (B4) 

for the remaining six 7r/2 rotations. For example, the action of the operators Rf on the Pauli operators/axes is 

Rt: iX,Y,Z)^iX,-Z,Y) (B5) 

R+: iX,Y,Z)^iZ,Y,-X) (B6) 

P+: {X,Y,Z)^{-Y,X,Z) (B7) 


Products of these nine elementary operations generate a representation of the 24 elements of the single-qubit Clifford group as 
tabulated in Table [In| We use this prescription to generate numerical simulations verifying our analytic calculations above. 
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# 

Gate Name 

Action on(X,Y, Z) 

Minimal Sequence(s) 

Notes 

Cl 

I 

iX,Y,Z) 

X^, R+R-, tG {1,2,3} 

identity 

C 2 

X 

{X,-Y,-Z) 

X 


Cs 

Y 

{-X,Y,-Z) 

Y 

TT rotation 

C4 

Z 

(-X,-Y,Z) 

Z 


C5 

Rt 

(X,-Z,Y) 

Rt 


Ce 

R+ 

(Z,Y,-X) 

i?+ 

-|-7r/2 rotations 

C7 

Rt 

(-Y,X,Z) 

R+ 


Cs 

Rx 

(X,Z,-Y) 

Rx 


C9 

Ry 

(-Z,Y,X) 

Ry 

—7r/2 rotations 

Cio 

R~ 

(Y,-X,Z) 

R- 


Cii 


(-X,-Z,-Y) 

ZR+ 


C12 


(-X,Z,Y) 

ZR~ 


Cl 3 


(-Y,-X,-Z) 

RtX 


Ci 4 


(Y,X,-Z) 

R-X 


Cl 5 


(-Y,-Z,X) 

RtR+ 


C16 


(-Y,Z,-X) 

R+R- 


Cl 7 


(-Z,-X,Y) 

R+R- 


Cis 


(-Z,-Y,-X) 

ZRy 


Cl 9 


(-Z,X,-Y) 

R+R- 


C20 


(Z,-X,-Y) 

R-R+ 


C2I 

H 

(Z,-Y,X) 

ZR+ 

Hadamard 

C22 


(Y,-Z,-X) 

R-Rt 


C23 


(Z,X,Y) 

R+R+ 


C24 


(Y,Z,X) 

R-RZ 



TABLE III. Representation of the Clifford group for a single qubit from products of elementary rotations. Relevant transformations of the 
coordinate system X,Y,Z under the action of each Clifford shown in column 3. Minimal sequence of elementary operations needed to 
generate each Clifford shown in column 4. We also indicate how these geometric rotations map to logical operations of interest for quantum 
information where relevant. 


C. Approximating the Noise-Averaged Fidelity {T) 


Here we present the full derivation of the approximation to the noise-averaged fidelity U^) given in Eq. 23 of the main text. 


Let 


denote the O 


(n 


xkp 

P=1 


= O (tr") term in the expansion of Eq. 


due to cross-multiplying terms like 


, where Retaining only terms up to fourth order (n = 4), consistent with our 

Talyor approximation of the error unitaries, we thereby obtain 


S,..«('■«> +d'’ + +fi5,.+ (SI ++ (SI, + (Sip+(SI + (S'’ 


(Cl) 
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where 


= Cl,. 




i=i 


C = Y.^^5i){i5k)Ci,,-iZCj,k.iZCk,j 

j<k 

^2^^ =~IJ2 S]Cij-iZ'^Cj^J 


i=i 


dll = Y. {i5j){^Sk){i5i)Ci^j-iZC,,k-iZCk,i-iZCij 

j<k<l 

j<k 


?3 

i=i 


j<.k<.l<.m 

11,2 = -\ E (*^j)(*4)<5;'Ci,,_iiC,-fe_iZCfe,,_iZ2C',,j 

j<k<l 


- ^ E {iS,)5l{i5i)Ci,,.iZCj,k-iZ^Ck,i-iZCi,j 

j<k<l 

-I Y S]{iSk){i6i)Ci,,.,Z^Cj,k-iZCk,i-iZCi,j 


j<k<l 


11 = 


j<k 

i 




j<k 

J 




i=i 


and we have defined the Clifford subsequence operators 

l<J<k<J. 

To evaluate Eq. |^we must obtain an expression for iTr(5,j,5). By the linearity of the trace, we must therefore obtain expressions 


(C2) 

(C3) 

(C4) 

(C5) 

(C6) 

(C7) 

(C8) 

(C9) 

(CIO) 

(Cll) 

(C12) 

(C13) 

(C14) 

(C15) 

(C16) 


for the quantities QiE 2 .....fc„ = [llM,...,k. 

highlight the following useful properties 


for each of the terms above. This requires some manipulation. To begin we 


= 


I Z, k odd 
I, k even 


Cl,J = Sr, =1. 

In fact since Ci, j = I, any cyclic permutation of subsequences Cjk also gives the identity. For instance 

Cl j = C\ j—\Cj k—iCk,j — Ck,jCi j—iCj k—i — Cj k—iCk,jCi,j—i ~ 2 


(C17) 

(C18) 


(C19) 
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and so on for any number of subsequences. Using the cyclic property of the trace, we therefore have 


Tr 




=-Tr 






iTr(z‘) 


0 , k odd 
1 , k even 


(C20) 


The last equality uses Eq. C17 the property of Pauli matrices that Tt(Z) = 0, and that the corresponding (2 x 2) identity has 


trace Tr(I) = 2. Now define the cumulative operators K^, giving the product from the first through to the mth Clifford operator 
in the sequence 


Krn = Ci,m= Cr,, ...Cr,^ , Ko=l=Kj (C21) 

where each Km is some element of the Clifford group due to the closure property of groups under group composition. In this 
case any subsequence Cjk “factorizes” into the products 

^ Kk = cl^_ .Crj- _ 1 Crj^ ■ ■ -Crf^. = Cj k (C22) 


allowing us to rewrite 


Cij_iZCj^k-iZCk,j = PjPk 
Ci,j.iZCj,k-iZCk,i-iZCij = PjPkPi 
ClJ-lZCj^k-lZCk,l-lZCl^rn-lZCm,J = PjPfePiPn 


where we define the .^-conjugating operators 


Pm = Km-lZKl_^ e {±X,±Y,±Z}, 


0 <m < J 


(C23) 
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and where is always a signed Pauli matrix due to the property that the Clifford group is the normalizer of the Pauli group. 
Thus we find 


= 1 (C24) 

= 0 (C25) 

=-I (C26) 

j<k 

Q2^ = -li2^J (C27) 

i=i 

Qili = -l Sj6kSiTv{P,PkPi) (C28) 

j<k<l 

Q[^1 = 0 (C29) 

Q® = 0 (C30) 

Qlll,l = l Sj6k6l6raTv{P,PkPlPm) (C31) 

j <.k<.l<m 

Q[% = \T. (PjPfe) + (P,P,) + S^SkSiTi (P^PO] (C32) 

j<k<l 

Q2l = \T.^Pk (C33) 

j<k 

^T3 = ^E^^^fcTr(P,Pfc) (C34) 

j<k 

= (C35) 

i=i 


The nonzero terms may be recast into more convenient expressions by moving to vector notation. We expand the operators P^ 
in the basis of Pauli operators, 

Pm — “t” yvaP “t” t ^1717 1/m5 ^ {O5 \^m \ “t” |l/m | “t” \^m \ — 1 - (C 36 ) 

That is, only one nonzero coefficient Xm,y,Zm, equivalent to expressing the fact that they are sampled from the set 
{±X, ±F, zLZ}. The associated unit vector 

t’m = {Xra7 ym^ Zm^ 7 Ijt^mH — 1 (C37) 

therefore points uniformly at random along one of the principle Cartesian axes {±5;, ±y, ±z}, capturing the “direction” of the 
operator Pm in “Pauli space”. With these definitions we may derive 

iTr(P,Pfc) =f,-f-fc (C38) 

^Tr (P.PfcPi) = i (f, X f-fc) • ri (C39) 

iTr(PjPfeP;P„) = {Vj ■ Tk) (r/ • T-m) + {rj ■ Tm) {Vk ' Ti) - (Vj ■ f/) (ffc • f„) 


(C40) 
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following directly from the trace and the cyclic composition properties of the Pauli matrices ai&j = ieijk^k + 5ijl, where Cijk 
is the Levi-Civita symbol, 6ij is the Kronecker delta and Einstein summation notation used. Consequently we obtain 

= (C41) 

j<k 

(C42) 

Om,i = X! (vj X fk) ■ ri (C43) 

j<k<l 

Qia,i,i = ^ X! Sj6k6i6m{(.rj ■ rk) {ri ■ Vjn) + (vj ■ rm) {rk ■ ri) - {rj ■ ri) {rk ■ Vjn)} (C44) 

Qm,2 = ^ X! ■ rk +SjSlSiTj ■ri+S’^SkSiT-k-ri] (C45) 

j<k<l 

Q^il = ■ r-k (C46) 

j<k 

Qi'hlT.^Pk (C47) 

j<k 

= (C48) 

f=i 


(O') 

We expect the major contribution to the distribution of (J^) to reside in the term Q\ ( since, in the limit a 1, higher order 

terms will be orders of magnitude smaller. Anticipating this, we recast in terms of an unrestricted sum to facilitate 

more straightforward analysis. Observing the quantity SjSkVj ■ Vk is invariant under exchange of indexes we thus obtain 


= - 4 (EE 

\j=l k=l 

= ■ Vk + 


j J 


where we define the resultant vector 


Hence 


-Tr 

2 


i=i 


2f.. . f . 

j' 3 '3 


j,k 


52 

2^ ^ 
i=i 


= “ 9 1 E^j^j I ■ (I 


( 2 ) 

2 


\k^l 


= -^\\Rr-Q\' 


R = Y.^3^3- 

i=i 


( 4 ,,) « 1 - - ll^f + ^ q(4)^^ ^ q(4) ^ + qO 


( 4 ) 


(C49) 

(C50) 

(C51) 

(C52) 

(C53) 

(C54) 


Since the above terms are all real, the fidelity is then obtained taking the square of Eq. |C54 and truncating any cross terms higher 
than O (cr‘^)- Thus we obtain 


J-(T7,<5)«l-||,Rf + 


fill 


+ 1 , 1,1 T ^Qi,i,2 ^^1,3 2^2,2 2(5^ C 


l( 4 ) 


l( 4 ) 


( 4 ) 


l( 4 ) 


l( 4 ) 


(C55) 
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Averaging Eq. C55 over an ensemble of noise realizations S then yields the approximate expression for (J^). However we 
make simplifying observation that only terms raised to even powers, or those summed over terms raised to even powers, survive 
the ensemble average. Moreover, since the vectors are uniformly distributed over the set {±x,±y,zLz}, summing over 
compositions 


{fj • f-fc) j <k < J (C56) 

{vj X fk) ■ fi j < k < I < J (C57) 

(fj • ffc) (f/ • Tm) j < k < I < m < J (C58) 

will on average resolve to the zero vector. Hence, even in the presence of correlated noise random variables, it is appropriate to 
set 

(C59) 

(C60) 

( 0 II 2 ) ^ 0 (C61) 

^ 0 (C62) 


In fact these quantities are random variables sharply (and symmetrically) peaked around 0. However since they are O (cr'^) 
terms, the spread of their distributions is very small relative to the spread of (||.R|p), which scales as O (tr^). The noise- 
averaged fidelity to O (cr"*) therefore reduces to 


where the final term in Eq. 


C63 


(j-)«i-(||^ir) + o(4) 

consists only of O (u^) terms 

oW = ^ (lifif) + 2 (qW) + 2 . 


(C63) 


(C64) 


The remaining task is to study how correlations between the noise variables affect the distributions of the terms in Eq. |C63[ and 
hence the distribution of (T). In particular how these distributions change depending on whether the noise falls into Markovian 
or quasi-static regimes. As stated in the main text, our key insight is the interpretation of the vector quantity R in terms of a 
3D random walk generated by adding J randomly-oriented steps with step lengths specified by S. The distribution of {R) then 
maps onto the distance square of this 3D random walk. 
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D. Universal Error Model 


Here we sketch a generalization of our analytic framework, showing its applicability to universal errors, beyond the (perhaps 
most important) case of dephasing specifically treated above. As in the main text, errors are implemented by interleaving Sri 
with a sequence of stochastic unitary rotations, yielding the noise-affected operation 


Sn.s = UiCniU2Cji2...lUjCn,. 
However we now let the unitaries take the general form 

Uj = exp 


-iAj ■ (T 


= exp 




■A<f'>ay + A[^ 


= exp 


-i (^A^pX + 

{x) \{v) 


(Dl) 


(D2) 


where cr = [ax,Oy, Uz) = (-A, Y, Z) denotes a vector of Pauli matrices and the vector Aj = {Aj, A^^\ ^j) contains the 
error rotations induced about each Cartesian axis. That is, the unitary Uj causes the net rotation through an angle || Aj|| about 
the axis Sj = A_,/|| || on the Bloch sphere. In this case the error process indicates a list of 3-component vectors 

<5 = (Ai,A2,...,A;). (D3) 

Assuming the weak-noise limit, quantified by the perturbative condition JE[|| A|p] < 1, we make the Taylor approximation 

Cj ~ I + * + A^- -f ... (D4) 

In the main text we presented a complete treatment up to fourth order to account for both the leading- and higher-order contribu¬ 
tions to infidelity, and assuming a dephasing-only environment. Here we demonstrate rigorously how leading-order contributions 
to the fidelity take the same form when the noise model is universal, as described above. Once again, the leading order contribu¬ 


tion derives from the term like Eq. 10 but with a slight amendment; 


where 




- ^t^C^,3-l^^.C,,k-lijuCk,J■ 


(D5) 


(D6) 


Writing Cju = and using Kq = Kj = I, we therefore obtain 

~ — 1 — 1 —1^1/ — 1 

j k 

where the operators in the last equality generalize the definition in Eq[T^ namely 

= Km-iairKl_^ e {±X,±Y,±Z}. 

Once again, 0 < m < J, and the Pm^ are always signed Pauli operators due to the property that the Clifford group is the 
normalizer of the Pauli group. We may therefore expand the operators in the basis of Pauli operators by writing 

P^)=a;^)A + 2/^)f + 4^)i (Dll) 

where Xm\ym \ Zm'^ S {0, ±1} with only one nonzero coefficient. Since the Clifford sequences (and hence subsequences) are 
uniformly random, the operators Pm^ are uniformly random, independent of the choice of y. The associated unit vector 

-(m) = (Xu-) ,,(m) Ilf(t*)II _ 1 

* m — v'^m •) Um ’ /'> W’ m W ^ 


(D7) 

(D8) 

(D9) 

(DIO) 


(D12) 
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therefore still points uniformly at random along one of the principle Cartesian axes {±5;, ±y, zLz}. Taking the trace over Eq. 
|D7| using the cyclic composition properties of the Pauli matrices, and moving to vector notation we therefore obtain 



(D13) 


Consequently 



(D14) 


j<k 


~ / 2^ 

where leading order contribution derives from the sum of all such terms Q\ { (/x, u) 


3 

Qm= E (D15) 

Taking an ensemble average over the noise realizations, we obtain 

{Q?i) = t. {q?1m) (D 16 ) 

11,iy—1 

= - E E (Af (D17) 

j<k 

We now make a simplification in which we assume the errors in separate quadratures arises from distinct physical mechanisms 
{e.g. dephasing and depolarization noise arise from different mechanisms). In this case the correlation between separate error 
quadratures is zero, and 



(D18) 


where is the Kronecker delta. With this assumption, we obtain 


3 



(D19) 


fj,=ij<k 


Now the quantity • f is invariant under exchange of lower indices, and we may recast the ordered sum as an 

unrestricted sum 



(D20) 


Defining 



(D21) 


we obtain 



(D22) 
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where the residual term 




(D23) 


^—1k—1 


is independent of the Clifford sequence and cancels out in the full expansion of (T), as in the case of the dephasing error model 
in the main text. The leading order contribution to the noise-averaged infidelity is therefore given by 




11=1 




(D24) 


where each of the terms (lli? 


ii2\ 


are independent random variable inheriting the PDF of the random walk, and taking the 

same form presented in the main text. Consequently each of the (|j.R p) are mutually independent, gamma-distributed 
random variables. In general, each of these random walks will be distinct, and the sum in Eq. D24| is over non-identical gamma- 
distributed random variables. The total PDF may be obtained from successive applications of Eq. |A4| by direct integration. In 
the event that all noise quadratures follow the same gamma-distribution r(a, fi), the total PDF simply reduces to the rescaled 
gamma distribution r(3Q;, /3) by an application of Eq. A6 
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E. PDF Derivation - Markovian Regime 


In the Markovian regime, we assume all noise random variables are i.i.d. Hence R corresponds to a J-length unbiased random 
walk with step lengths sampled from the normal distribution JV (O, cr^). Since these step lengths are symmetrically distributed 
about zero, the distributions of the components of the walk vector SjTj = {SjXj, Sjyj, SjZj) are invariant with respect to the 
sign of the coefficients aj in all Cartesian directions a G {cc, y, z}. Ignoring the signs we therefore treat the the coefficients as 
binaries aj G {0,1}, where the zero event simply reduces the number of steps taken in that direction. Let 

.7 

na='^\aj\, aG{x,y,z}, ria: + riy + Uz = J (El) 

count the total number of nonzero components in each Cartesian direction over the sequence of walk vectors {fi, f 2 ,f j}. 
Thus 


+ + sy + ,„ + sy^, + + (E2) 

where the superscripts in (5“ indicate summing only over the subset of 6j for which the coefficients aj are nonzero. Thus we 
have 


\\Rf = Al + Al + Al A„= + ()“ + ...+ J“J, aG{x,y,z}. (E3) 

Since all 6j ~ Af (O, a^) are i.i.d. in the Markovian regime, so too are all the random variables in set Sa = ^2 j 

a G {x, y, z}. The distribution of their sum is therefore given by 

- AA (0, ric^a^) . (E4) 


Eurther, since each f-j projects onto only a single Cartesian direction - and consequently the sets Sx,y,z are mutually disjoint 
- the random variables A^^y^z are mutually independent. The distribution of the sum of squares of non-identical Gaussians 
involves a generalized generally challenging to write down, requiring a generalized chi-square distribution. To preserve analytic 
tractability, we make the following simplihcation. Since the vectors rj are uniformly-distributed there is a | probability of being 
parallel to any given Cartesian axis. The probability of getting any particular combination {nx, Uy, riz) is therefore given by the 
multinomial distribution 




1 ^ 7 / ? 


:) = 


J}_ 

Ux'-riylnz'. 




Ux + riy + riz = J 


(E5) 


which, for J ^ 5, is sufficiently peaked around nx,y,z = J/^ that we may regard these values as fixed without signihcant 
error. In this case Ax^y,z ~ Af (O, Jcr^/3) reduce to i.i.d. random variables. The distribution of ||.R|p consequently reduces 
to chi-square distribution with 3 degrees of freedom. It is more convenient, however, to express this in more general terms as a 
member of the two-parameter family of gamma distributions (see Eq. |A5| l, of which the chi-square is a special case. Specifically, 
we obtain the gamma distribution 


l-RII 


r {a, 13), 


2Ja^ 


“=2’ 


(E6) 


with shape parameter a and scale parameter /?. The distribution of a hnite noise-ensemble average over ||i?|p is therefore 
specified by 


Rift.. = rtf iiRiij-rfs."^) 


(E7) 


where now the || ii||^ are i.i.d. random variables. But the sample mean over n gamma-distributed random variables simply yields 
a rescaled gamma distribution with a -G na and [3 ^ (3/n (see Eq. A7 1 . Consequently 

mw 


"3n 2Ja^ 
T’ 3n 


(E8) 
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with moments 


E[(|l^f )„] = Ja2 

miiRfu = 


(E9) 

(ElO) 


From Eq. 


E6 


the PDF for liijlP now has the known form 




T{a)P' 


-a;“ ^ exp 


(Ell) 


where a = 3/2, /3 = 2 Jct^/ 3, and r(a;)is the gamma function. The PDF for ||ii||^ therefore given by the transformation (see 
Eq.[A^ 


1 


(E12) 


By direct computation, the first two moments are then given by 


E[\\Rf] = - JV 
V[|1.R||^1 = jV 


(E13) 

(E14) 


It is relevant at this point to consider the relative weight of the terms (|ji2|j^)„ and (|jil|p)„ in the calculation of (T). An 
application of the central limit theorem yields the approximation 


■ ((||fir)„ - E[(||^r)„]) ^ U (o, Y[\\Rf]) 


That is, (||i?||‘*)„ is approximately Gaussian-distributed with mean and variance 


(E15) 


E[{\\R\\%]=E[\\Rf] = -J^a^ 

RO 

V[(|l^r)„] = V[||,Rr]/n = - jVn-i 


(E16) 

(E17) 


The relative significance of these terms may be captured by the condition V[(||i2|| ^)n] < eV[(||.R|l2)„] where e is some small 
fraction. From Eqs. ElO and E17 this condition is met provided 

J< 


(E18) 


-2 


For instance, if cr ~ 0.01, e < 0.1 provided we restrict J < 1000. Now, since {R) involves a linear combination of the 
terms (||ii||'‘)„ and (||i2|p)„, the PDF of (R) is approximately given by a convolution over the individual PDFs of these terms. 
However assuming e is sufficiently small - so that the distribution of (||.R||"‘)„ is sufficiently narrow - the primary contribution 
of this convolution is to shift the distribution by an amount approximately given by the mean of (||.R||^)„. Thus we set 


(ii^r)„^E[(ii,Rr)„] = -jv 


(E19) 
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The other O (a^) terms reduce to 


\j<k / j<k 




\i=i 


i=i 


(E20) 

(E21) 


Substituting in these constants, Eqs. |C63|and|C64| the noise average fidelity for Markovian errors reduces to 




{\\Rr)n-T 


3n 2J<7^ 


2 ’ 3n 


(E22) 


Performing the appropriate linear transformations (Eq. |A1[ ) to incorporate the constant factors in Eq |E22| and using the definition 
of the PDE for a gamma-distribution (Eq. A5 1 , the PDE for {R) finally takes the form 


i'(r) = 1 - F + 

3 


.(F) 


= 3n/2, /3 = 2Ja^l3'r 


(E23) 

(E24) 

(E25) 


where r(a;)is the gamma function. 
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F. PDF Derivation - Quasi-Static (DC) Regime 


In the DC regime we assume all noise random variables 6j = 6 are, in a given instance, identical (maximally correlated). 
However over separate instances S is sampled from the normal distribution 6 ^ Af (O, In a given run, the random walk 

vector R — ^ therefore reduces to a J-length unbiased random walk on a 3D lattice with fixed step length 

S. In this case the noise random variables S and Clifford-dependent random variables rj factorize, allowing us to express 


R = SV, 




(Fl) 


i=i 


where V G defines an unbiased random walk on a 3D lattice generated by adding J unit-length steps. Since we are interested 
in the norm square ||i2|p = i5^|| V|p, however, any sign dependence of 5 vanishes. Performing a finite ensemble average over 
noise randomizations we therefore obtain 


\Rf)n = {SX\m 


{\\R\\^)n = {S^)n\\V\\ 


(F2) 


Substituting these expressions into Eq. C63 we therefore obtain {T) « 1 — (<)^)n||^|P + ^ where the 

term O (cr^) includes the quantities (Q 2 2 ) (^ 2 ^ 2 )’ which approximately reduce to constants. Since the term in || and 

IIV11^ are now highly correlated, however, we cannot exploit simplifying properties such as approximate independence between 
primary and higher-order terms in the expansion. One approach would involve making the approximation that {5‘^)n and 

for sufficiently large n, and then completing the square in || V|l^. In this case the PDF for {R) could be obtained 
by successively performing linear, square and shifting transformations on the PDF of || Vp. However to a good approximation 
it turns out most of the physics is captured by the first term (i5^)„|| Vp. Hence we make the approximation 


(J-)«1-(<52)„||F|| 


(F3) 


As shown in the main text, proceeding with this truncating the expansion produces good agreement with direct simulation. The 
PDF for {JF) is therefore obtained by incorporating the distributions of || V|p and T he PD F of 5 is given by the Gaussian 

fs{x) := exp ■ Hence the PDF for is given by the transformation (see Eq. 


A31 


fs 2 (x) = 2 





exp 




r(i)(2a2)i/2' 


r(a)^‘ 


^ exp 


exp 

X 

1 




2a^J 


(E4) 

(F5) 

(F6) 

(F7) 


where a = 1/2 and j3 = 2a^, and r(a;) is the gamma function. Comparing this with Eq. A5 this is the PDE of a Gamma 


distribution with shape parameter a and scale parameter /3. The distribution of a hnite noise-ensemble average over 62 is 
therefore specified by 




1=1 


r( 2 . 


(F8) 


where the Sj are i.i.d. random variables. Now the sample mean over n gamma-distributed random variables simply yields a 
rescaled gamma distribution with a —>■ na and ^ fijn (see Eq. A7 1 . Consequently 




I'i-- 


(F9) 
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with moments 


E[((52)„] = 

V[(<5^)„] = — 


(FIO) 

(Fll) 


As outlined above, ||V|p expresses the distance square of an unbiased random walk on a 3D lattice generated by adding J 
unit-length steps. Let R be the random variable representing the distance from the origin in a symmetric (Bernoulli) 3D random 
walk after J steps. Then the PDF for R is known to be 


= (2I7) 


3/2 


9 

47rr e 


The distribution of the distance square || V|p = R^ is therefore given by the transformation (see Eq. 


A3 I 




2x1/2 


/« 


3/2 

-3/2 

-3/2 


= %7rxi/2 

\27rJ J 


2J 

T 

2J 


exp 

-1 


—3x 


= ^ r ^ 


„a-l 


x^/^ exp 


x^^^ exp 


— X 

2J73_ 

— X 

2J/3 


r(a)/3° 


expl-^ 


(F12) 

(F13) 

(F14) 

(F15) 

(F16) 

(F17) 


where a = 3/2 and f3 = 2J/3, and r(x) is the gamma function. But this is the PDF of a Gamma distribution with shape 
parameter a and scale parameter /3. Consequently 


with moments 


IVIi 


Fl^ ^ 

2’ 3 


E[\\Vf]=J 

V[||Vf] = ^ 


(F18) 

(F19) 

(F20) 


Thus, to first order the PDF for {T) is specified by the product of independent gamma-distributed random variables. This class 
of distributions is generally difficult to express. However we may obtain a closed form for the PDF of (^^)n||V|p by direct 
integration (see Eq. |A2|i, yielding 


f{r) - 


n + 3 
K 4 



Ki 


(n-3) 



(v/kF) 


(E21) 


where v = 1 — F, n = Sn/Ja"^ and Kn{z) gives the modified Bessel function of the second kind. However for fairly reasonable 
ensemble sizes n > 50 it is sufficient to approximate as the mean of its distribution, namely cr^. In this case the fidelity 
distribution reduces simply to the (scaled and shifted) gamma distribution associated with which || V|p, yielding 


v=l-F 
a = 3/2, 


^ exp 



^ = 2JcrV3 


(E22) 

(E23) 

(E24) 
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G. Probability distribution functions for partially correlated noise 

In the main text we derived gamma distributions for (J^) in the Markovian and DC limits, and assumed continuity of the 
distribution to interpolate between these extremal cases. In this appendix we present formal justification of this approach. We 
treat a specific intermediate-correlation-length model where errors are block-correlated, i.e. within blocks of length M, the noise 
value is identical, and there is no correlation between blocks. That is, each noise realization, 6, can be partitioned into blocks 
of length M, within which the Si ^ JV{0, a^) are constant. In block k, the corresponding random walk takes steps along the 
cartesian axes (in either direction), which we count as counts satisfy the constraint 


“ -I- + Wfc ^ = M 


and are multinomially distributed, 

— 1^+^ ryy^-^ 




(Gl) 


(G2) 


The displacement vector associated with block k is then 


( 


- rrit. 


= SkVk 


(G3) 


where Vk is. the displacement vector of a random walk associated with the fcth block, involving M unit-length steps along the 
cartesian axes. The total displacement is then 


N 

R = '^Rk, N = J/M (G4) 

/c=l 

where N is the number of M-length blocks in a sequence of total length J. In this picture the Markovian limit corresponds to the 
case M — 1,N = J {J blocks, each consisting of a single Clifford gate); the DC limit corresponds to the case M = J, N = 1 
(1 block of length J). We therefore obtain 


N N 


= EE^*-^^- 

N N 

(G5) 

= Eii^fef 

k^l i^j 

N N 

(G6) 

= j2^k\\Vkr+j2s.s,v.-v,. 

(G7) 


k=l i=fj 


Taking the ensemble average over noise realizations then yields 

AT N 

(II.RI12) = J2{Sl)\\VkV + 


(G8) 


Here the angle brackets (•) refer to an ensemble average over noise realizations in the limit of an infinite ensemble size (n —> c»). 
From our assumption that the noise is perfectly block correlated {i.e. zero correlation between errors from different blocks), and 
Gaussian-distributed within each block Si ^ JV{0, ct^), we therefore obtain 


(SiSj) = aHi 


(G9) 


where Sij is the Kronecker delta. Thus, the second term in Eq. G8 vanishes, yielding 


N 


N 


(ii^f)=E E - mrr + - mr?)] 


(GIO) 
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In general, (||i2|p), depends on the details of the Clifford sequence rj through its dependence on the Vk, or equivalently, on 
the step counts mt. Considering an ensemble of Clifford sequences, we can compute the statistics of (||i?|p). The first two 
moments are straightforward to calculate exactly. We find 


E[((m+^ - + (m+^ - + (m+^ - = a'^M, 

V[(mr - = 2 a^M{M - l)/3, 


which, together with Eq. |G10 therefore yields 


EKIlHf)] = aV, 
V[(||.Rf)]=2aV(M-l)/3. 


(Gil) 

(G12) 


In the Markovian limit (M = 1,7V = J), the expectation E[(||ii|p)] = a^J is in fact independent of the details of the 


sequence, since m consists of a single unit entry, so that the summand in Eq. |G10 is identically 1. This is perhaps not surprising: 
if different Clifford sequences had different noise-averaged fidelities for Markovian noise, then the optimal sequences would be 
good candidates for dynamical-decoupling sequences to suppress Markovian noise. It is known that dynamical decoupling is not 
useful for correcting fast noise. The variance is consistent with the n —>■ oo limit, V[(||i?|p)oo] = 0 in the main text (see Table 
[I]l. Eor DC noise (M = J,N = 1), we find the same Clifford-averaged expectation as above, however 


)] = 2(tV(J - 1)/3 « 2(7^JV3 


(G13) 


in agreement with earlier results in this limit up to a term O (cr^) (see Table^. 

In a given experimental scenario, it may be that there is a maximum corrmation length Mmax, beyond which gates are uncor¬ 
related. The above results then suggests that the clifford-averaged variance will vary quadratically with J for J ^ Mmax, and 
linearly with J for J ^ with a transition around J = Mmax- This may be a useful heuristic for identifying correlation 

lengths in experiments. 

For a given sequence length, J, the PDF of (||.R|p) is in fact discrete. However, for modest sizes of J and M we can 
approximate the PDF by a continuous distribution. Given that the PDFs in both the DC and Markovian limits are both well- 
approximated by gamma distributions, we can use the calculated expectation and variance above t o gues s t he PD F in the regime 
of intermediate correlation length. Since E[r(a, /3)] = aP and V[r(a, /?)] = a/3^, and using Eqs. Gil and G12 we guess 


|H||2)^r( 


3J 2(M - l)cr2 


2(M-1)’ 3 

/ 3J 2Ma^ 
\2M' 3 


for M > 1 


(G14) 

(G15) 


Eq. G15| interpolates between the two limiting cases. 

This guess can formalized by recalling (||.R|p) = However, each of the Vk represent the displacement of 

a independent random walk taking M steps of unit length along the cartesian axes, and are formally identical to the random 


walk vector in the DC limit. From Eq. 41 the displa cemen t square of this random walk, given J steps, is given by the gamma 
distribution |j ~ T (|, ^). Consequently Eq. GIO describes the sum, scaled by of N independent and identically 


distributed random variables, each following the distribution 

llV.f- 

By an application of Eqs. |A6| and [A7| we therefore recover 

m?) 


2’ 3 


\2M' 


2Mcr" 


(G16) 


(G17) 


which is asymptotically correct for large M. In showing this we have demonstrated analytically that intermediate correlation 
structures between the Markovian and DC limits are also described by the gamma distribution. 



















30 


H. Fidelity Statistics for Generically-Correlated Processes 


Here we derive the expectation E[(||i?|j^)] and variance V[(||fi|p)] stated in Eq. 45 of the main text, for an error process S 


with generic correlation structure specified by an autocorrelation function. Consider an arbitrary time-series x{t) describing 
some (continuous) time-dependent, wide-sense stationary error process. Then x{t) may be characterized by an autocorrelation 
function 


C^{t) = {x{t)x{t + r))j = ^1^ ^ 

where ()t refers to an ensemble average over the time series and r is the time difference between measurements. Invoking the 
Wiener-Khintchine theorem, Cx{t) and the power spectral-density (PSD) S{uj) form a Fourier-transform pair 

-| noo poo 

S{uj) = / C'(r)e“^dT C(r) = / S{uj)e-^^^duj (H2) 

We make use of these relations by discretizing the time-series, whereby 

j 

R = '^djrj, 6j=x{tj), tj/Tg e {1,2,.., J} (H3) 

and Tg is the time taken to perform a Clifford operation. The underlying error process is thereby discretely “sampled” by the 
Clifford sequence Sr/, and correlations between elements of d separated by a time interval of “k gates” are specified by the 
(discrete) autocorrelation function 


Cs{k) — {SjSj+k)j 


(H4) 


Expectation 


We begin by expanding 


,7 ,7 

(IIAf) = {R-R)= ■ r,. (H5) 

i=i j=i 

Here the ensemble average indicated by (•) is over the noise random variables, and so does not affect the Clifford-dependent 
unit vectors f^. To obtain the full expectation we now take a separate expectation over the the random variables f^, which does 
not affect the noise random variables. Since the Clifford unit vectors f ^ are uniformly and randomly distributed over the set 
{±x, ±y, zLz} the expectation over the inner product E[fi • Vj] = Sij (Kronecker delta). Hence 


.7 .7 


E[(||-Rf )] = y]] • Vj] 

(H6) 

i=i 


.7 .7 


= 'y (y '{x(ti)x(tj))6ij 

(H7) 

i=l j = l 


J 


i=l 

(H8) 

.7 

(H9) 

.7 

= E^‘5(0) 

2—1 

(HIO) 

= JC5(0) 

(Hll) 
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Variance 


For the variance we need to find the expectation of (||ii|p)^. Using the above result, we momentarily define 

J ,/ 


£:^E[(|l,Rf)] = ^(44) = E(^fc) 


(H12) 






Now we rewrite 


\Rf) = ■ rj) = £ + ■ r,) 


(H13) 


k=l 


Hence 


l-Rf ))^ = ( £ + • fj) I I £ + ^ (4<5/)(ri/ • fjv) 

\ / \ i'/i' 

J J J 

= £^ + 2^2 ^((5,4)(fi • fj) + {SiSj){SeSf){r, ■ fj){re ■ r^,) 

i ¥=3 i ¥^3 


(H14) 


(H15) 


As above, E[ri • rj] = 6ij, so taking the expectation over the the random variables f ^ causes the term 
vanishes, yielding 


E[((|!,Rf ))2] =£^ + ^ fl { 5 .S,){S, 6 ,,)E[{r^ ■ r,)(r.' • Vf)] 

i'¥=3' 

Hence the variance is given by 

V[(|l,Rf)]=E[((|l,Rf))2]_E[((|lHf))P 

= E[((|l,Rf))Vf" 

J J 

= H • rj){r,, ■ Tj,)] 

i¥^3 i'/i' 

Now for i ^ j and i' ^ j', the expectation over (f^ • rj){rii ■ Tji) is zero unless either 

i = i' k, j = j' 
i = j' kj = i! 

The variance therefore reduces to 

mm")] = 

i¥=3 iVj *7^i 

Observing the correlations 

(Sidj) = {6iSi+(^j_,)) = Cs{j - i) 


E[(o •T'i)(0' -rf)] = 1/3 
E[(o •T'i)(0' -rj')] = 1/3 


(H16) 


(H17) 

(H18) 

(H19) 


(H20) 

(H21) 


(H22) 


(H23) 


depend only on the difference index k = j — i, which runs between —{J — 1) and J — 1, we may reexpress the sum over all 
i^j as 


2 

mm\")] = 3 E -1^1) - -^(^^(0))' 


(H24) 
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Since Cs{k) = Cs{—k), we may reexpress the above simply as a one sided sum 

^ k=l 

where we have removed the contribution J(C' 5 ( 0 ))^ by commencing the sum at fc = 1. 


(H25) 
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I. Generating Aribtrary PSDs from Fourier Synthesis 

Following previous work on engineering noise processes ||29|, we construct the x{t) as a superposition of phase-randomized 
cosines 


x{t) 


auto 


Q 

9=1 




iliq ~ UniformDistribution[0, 2tt] 


(II) 


Here a is a global scaling factor for setting the total power content, ujq is the mode separation incrementing the (angular) 
frequency in the Fourier superposition, and the mode frequencies are given by ujg = ui^q. The function F{q) specifies the 
relative weighting of the Fourier components. The autocorrelation function is then given by 

C'x(r) = {x{t + T)x{t))t = '^{qF{q)f cos ( w , t ) (12) 

9=1 


where (•)( denotes averaging over all times t from which the relative lag of duration r is dehned. Invoking the Wiener-Khintchine 
theorem and moving to the Fourier domain we then obtain the power spectral density 


S(uj) = '^{qF{q)f S(UJ - UJq) + 5{uj + UJq) 

9=1 


(13) 


Thus in this formulation the power-spectral density is represented as a Dirac comb of discrete frequency components with the 
amplitude of the jth tooth determined by the quantity {q{F{q)Y. It is then straightforeward to specify the construction of any 
power-law PSD by writing the amplitude of the qth frequency component as a power-law, S{oj) cx {qui^Y- It therefore follows 
that the envelope function for the comb teeth in the phase modulation scales as 


F{q) = q 2 


£_1 


Table IV shows the functional form required for F {q) in order to achieve dephasing-noise PSDs of interest. 


(14) 



Power Laws 

1//^ 1// White Ohmic 

p 

F{q) 

-2-10 1 
q-2 g-3/2 q-i q-i/2 


TABLE IV. Functional form of F{q) for well-known power-law PSDs 


From the standpoint of the Clifford sequence sampling the noise process, the PSD reconstruction S{ijj) is obtained as the band- 
limited Fourier transform of|^ restricted to the domain r G [—J, J]. Hence 


S{oj) = [0(r + J) - 0(r - J)] } 


^ 7 ^ [sinc( J(a; - ujj)) + sinc(J(w -f w^-)) 

' 0=1 


(15) 

(16) 


where 


lim S{uj) = S{uj) 

J —¥00 


( 17 ) 
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J. Confidence Bounds for 

In the main text we have showed the noise averaged fidelity {F) is, to first order, a random variable specified by 

{F) = A-iy (Jl) 

where ly ^ T (a, /3) is a gamma-distributed random variable, and the values of the parameters A, a and P are dependent on 
character of the noise. In particular, restricting attention to the Markovian or DC regime, we have 



Markovian 

DC 

a 

|n 

3 

2 

/? 



A 


1 


TABLE V. Input parameters for calculating PDF of (F). 

The true mean fidelity formally obtained as an expectation over all possible fidelity outcomes F defined on the support of the 
random variables t] and d is given by 


{F}r,^s = E[(.F)] = A - a/3 (J2) 

These equivalent expressions for the total expectation value are hereafter denotes simply by as in the main text. In the 

standard RB procedure this expectation value is estimated by the sample mean 



1=1 


{■) 


(J3) 


The mean gate error prb is then approximated by the decay constant from an exponential fit for increasing J. It is therefore of 
interest to quantify the reliability of the estimate as a function of the ensemble size. This may be expressed in terms of 

confidence intervals, treating the measured values as random variables sampled from the distribution describing {F). The 
distribution of is therefore specified by 




(■') = A - 17, 


k 


~ r (a,/3) 




(J4) 


where, since the Vi are i.i.d. random variables we have j/ ~ T {ka, P/k). The probability that falls outside the confidence 
interval j — L, E[J^] j + U] is therefore given by 


S=l-V (E[J'] J -L<A-V< E[J'] j -f U ), 


k 


(J5) 


where L,U > 0 specify the lengths of the lower and upper error bars centered on the expectation value E[J^]j. Substituting 
E[J^] J = A — a/3 and directly integrating the PDF over these confidence bounds we obtain 


6{a, P, L, U, k) 


l-Q 


ka, k 



,k 



(J6) 


where Q{a, zq, zi) is the generalized regularized incomplete gamma function. This is defined in non-singular cases by 


Q{a,zo,zi) 


r{a,zo,Zi) 

r(a) 


(J7) 


where r{a, Zq, zi) = r(a, zq) — r(a, Zi) is the generalized incomplete gamma function, r(a, z) is the incomplete gamma 
function, and r(a) is the Euler gamma function. This expression can be further condensed by scaling the error-bar lengths by 


L = GL{A-E[F]j) = GLaP 
U = Gu {A-E[F]j) = GuaP 


(J8) 

(J9) 
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where the values Gl, Gu > 0 specify the lengths of the error bars as fractions of the expected infidelity. In this framework, the 
reliability of the estimate may be quantified by the requirement that 

S(a, (3, GlQ-P, Guctfi, k) < e (JIO) 

where e is a small fraction. Substituting values of a and (3 into Eq. |J6] appropriate for Markovian (M) and DC regimes, we find 

= l-Q 


3kn 3kn ^ , 3kn. ^ , 

^(1-Gj/), —(1 + Gl) 


SiDC) ^ 1 _ Q 


3k 3k ^ \ 

—, -a-Gu), y(1 + Gl) 


(Jll) 

(J12) 


For user-defined Gl, Gu and e, one may then solve the inequality in Eq. |J10| for minimum k. Thus one may bound from below 
the size of the ensemble k necessary to justify the uncertainties quoted for the mean gate errors prb obtained from RB. 







